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Abstract 

The solid-liquid phase-diagram of charge-stabilized colloidal suspensions is calculated 
using a technique that combines a continuous Poisson-Boltzmann description for the mi- 
croscopic electrolyte ions with a molecular-dynamics simulation for the macroionic colloidal 
spheres. While correlations between the microions are neglected in this approach, many-body 
interactions between the colloids are fully included. The solid-liquid transition is determined 
at a high colloid volume fraction where many-body interactions are expected to be strong. 
With a view to the Derjaguin-Landau-Verwey-Overbeek theory predicting that colloids in- 
teract via Yukawa pair-potentials, we compare our results with the phase diagram of a simple 
Yukawa liquid. Good agreement is found at high salt conditions, while at low ionic strength 
considerable deviations are observed. By calculating effective colloid-colloid pair-interactions 
it is demonstrated that these differences are due to many-body interactions. We suggest a 
density-dependent pair-potential in the form of a truncated Yukawa potential, and show 
that it offers a considerably improved description of the solid-liquid phase-behavior of con- 
centrated colloidal suspensions. 



The classic Derjaguin-Landau-Verwey-Overbeek (DLVO) theory || predicts that an isolated 
pair of charged colloidal spheres in an aqueous salt solution interacts via a repulsive Yukawa 
potential at large separations. Direct experimental measurements have confirmed the general 
validity of this prediction pL ||, [| . Though originally being designed only for pairs of colloids in 
isolation, it is common practice nowadays to apply the DLVO theory also to the much more com- 
plex case of colloids in concentrated suspensions, the situation which occurs in most technical 
applications. As we show in this letter, treating a colloidal suspension as a system of pairwise 
interacting Yukawa particles is valid only in a limited parameter regime. To this end a full 
Poisson-Boltzmann (PB) mean-field description of the colloidal interactions has been combined 
with a molecular-dynamics (MD) simulation, similar to Ref. ||, to calculate the solid-liquid 
phase-behavior of charge-stabilized colloidal suspensions. While correlations between the mi- 
croions are neglected in this approach, many-body interactions between the colloids, mediated 
by the screening ionic fluid between and around the colloids, are fully included. By calculating 
the effective DLVO parameters for the same system, it is seen that at large volume fraction 
of colloids and low salt concentration, the description in terms of pair-wise interacting Yukawa 
particles fails dramatically. 
1 Jure.Dobnikar@uni-konstanz.de 
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An estimate whether the DLVO pair potential is a good approximation in concentrated 
colloidal suspensions can be based on the ratio between i) the mean colloid-colloid distance, 
<^m = p 1 / 3 ) in suspensions with a density of colloidal particles p and ii) the inverse Debye 
screening length k^ 1 , which measures the thickness of the spherical double-layer around a single 
isolated colloidal sphere due to screening by the microions. If two or more such spherical 
double-layers overlap, the screening of the colloidal charges becomes incomplete and the charge 
distributions on the colloids involved begin to interact. If the colloid density is low, i.e. if 
d m K 3> 1, it is obvious that such an overlap will almost always occur for pairs of colloids only. 
At high colloid densities, however, i.e. if d m K ~ 1, there is a high probability that more than 
one other colloid is within the range of the double-layer k _1 around any colloidal particle. In 
this case many-body forces between the colloids become important, and the picture of pair-wise 
interacting particles ceases to be valid. 

This scenario has been qualitatively confirmed by recent measurements of the effective col- 
loidal pair-interaction potential as a function of the colloid density in 2D colloidal systems 
under low-salt conditions j|, g. In this context, many-body interactions show up as a density- 
dependence of the effective pair-potential 0. Indeed, at low density, the measured pair-potential 
was Yukawa-like as predicted by DLVO theory, but at high density it showed clear deviations 
from a Yukawa-like behavior, particularly for inter-colloidal distances r > d m where the in- 
teraction decayed to zero much more rapidly than predicted by DLVO theory. This observed 
density-dependent "cut-off" of the interaction potential at the mean colloid-colloid distance 
d m is a direct manifestation of many-body interactions in the system ||. A simple physical 
explanation for this "cut-off" is based on the mechanism of macro- ion shielding ||: the in- 
teraction between two colloids separated by a distance r > d m is likely to be screened by a 
third colloid located somewhere between them. This screening effect of the macroionic charges, 
which is neglected in the DLVO theory, qualitatively explains the observed deviations from the 
Yukawa-potential. Motivated by these experiments, we here address the question to what extent 
many-body interactions affect macroscopic properties of 3D colloidal suspensions, specifically, 
the solid- liquid phase behavior. 

Colloidal crystals near melting are simulated using a technique suggested by Fushiki [g] which 
combines a PB field description for the microions with a MD simulation for the macroionic col- 
loids Q. Specifically, we consider a fluid of highly charged identical colloidal spheres suspended 
in a structureless medium of dielectric constant e at temperature T. The Bjerrum length char- 
acterizing the medium is defined as Xb = e 2 P/e, with (3 = 1/kT, k Boltzmann's constant, and 
e the elementary charge. The colloidal spheres have radius a, charge — Ze, and thus a surface 
charge density —ea = — Ze/Ana 2 . They are placed in a cubic box with periodic boundaries and 
the positions of their centers are denoted by Rj (i = 1 . . . N). The suspension is assumed to be 
in osmotic equilibrium with an electroneutral reservoir of monovalent point-like salt ions with 
total particle density 2c s and inverse Debye screening length k = (8-kXbCs) 1 ^ 2 , and in thermal 
equilibrium with a heat bath at constant temperature. The density distribution of the elec- 
trolyte ions in the region G exterior to the colloidal spheres is obtained from the normalized 
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electrostatic potential <^>(r) satisfying the PB equation 

V 2 4>(y) = K 2 sinh0(r), r € G 

mV(f> = 47r\ B a r EdGi, i = l,...,N, 

where dGi is the surface of the i-th colloid with outward pointing surface normal rij. Constant- 
charge boundary conditions are assumed for all ./V colloid surfaces [^Oj] . The PB-MD simulation 
algorithm is now described as follows ||: (i) fix all N colloids at their initial positions {Rj}°, (ii) 
solve the PB eq. (|l|) for this colloidal configuration, (iii) calculate the total force on each colloid 
by integrating the stress-tensor over a surface enclosing the respective particle, (iv) move all 
colloids to their new positions by an Euler MD time-step including a stochastic force modeling 
the heat bath. Repeat (ii - iv). The solution of the PB equation in step (ii) is computationally 
the most demanding part of the calculation since it is necessary to resolve the steep variation 
of <p near the colloid surfaces. To achieve this resolution with a reasonable number of grid- 
points, N spherical coordinate systems are constructed, one centered around each colloid, which 
overlap with a Cartesian system in the cubic simulation box. The solution is obtained by (iia) 
solving the PB equation in each of the spherical systems for fixed values at the outer edges, (iib) 
interpolating to find the boundary values in the Cartesian system, (he) solving the PB equation 
in the Cartesian system, and (iid) interpolating back to find new boundary values at the edges 
of the spherical systems. These steps are iterated to convergence. There are a few studies 



TIL 12, 13, ra] where similar simulations have been performed, but they mainly concentrate on 



determining the structure of colloidal dispersions. The method used by Lowen et al. [11, [T^] is 
based on a density functional approach and includes also the microion correlations. However, it 



has been shown in various studies [14, 15 1 that the microion correlations play only a negligible 
role in case of monovalent salt ions as considered in this paper. Groot []l4| quantified the validity 
of the PB approach by comparing it to a cell model Monte-Carlo simulation and showed that 
the deviations are already tiny for a ratio Ab/cl = 0.03. Since we are here considering an even 
smaller value Xb/cl = 0.012, the PB approach is perfectly justified relying on Groot's work. 

Using the simulation method just described, we have calculated the melting line in the phase- 
diagram of charge-stabilized colloidal suspensions. To facilitate later comparison of our results 
with the classic phase-diagram of Yukawa systems calculated by Robbins, Kremer and Grest 
(RKG) [[IT]], we followed their work and determined the phase-boundary by the Lindemann rule 
which states that a crystal begins to melt when the rms displacement in the solid phase is a 
fraction of 19 % of d m . Taking the radius a of the colloidal spheres as unit-length scale and 
having in mind an aqueous system at room temperature where Xb = 7.2 A and colloid particles 
of size a = 60nm as common in experiments we fix the ratio Xb/cl = 0.012. The state space of 
our colloidal system then is three-dimensional, spanned by na (salt concentration), Z (colloidal 
charge) and d m /a (colloid density), or, equivalently, the volume-fraction r/ = &ir(a/d m ) 3 /3. A 
point X in state space is thus given by X = (Ka,Z,rj). Since a systematic exploration of the 
full 3D phase-diagram is too expensive computationally, we focus on a 2D cut at large volume 
fractions r/ = 0.03 realized in experiments. The resulting phase-diagram in the (Z, «a)-plane 
shown in Fig. (|l|.a) represents our main result. It describes the solid-liquid phase behavior 



3 




Figure 1: (a) Solid-liquid phase diagram of a charge-stabilized colloidal suspension, spanned 
by the colloidal charge Z and the salt concentration na for volume fraction r/ = 0.03 and 
As/a = 0.012. The circles are points on the melting line obtained by applying the Lindemann 
criterium to PB-MD simulation data in a FCC crystal which fully include many-body inter- 
actions between the colloids. The line is a guide to the eyes, (b) The melting line of a pure 



Yukawa system obtained by Robbins et al. [17] (solid line), the melting line obtained assuming 



Yukawa interactions with a density-dependent cut-off after the first neighbor shell in a FCC 
configuration (dotted line) and in a BCC configuration (dot-dashed line), and the data from 
Fig. (jjja) transfered by using the effective force parameters as described in the text. 

of colloidal suspensions with all many-body interactions among the colloids included. To our 
knowledge, it is the first such phase diagram calculated without assuming pairwise additive 
effective interactions. Its significance is best appreciated by contrasting it to the predictions 
based on the established description in terms of pairwise Yukawa potentials, which in turn 
requires to first introduce the concepts of charge renormalization and effective forces. 

For a system of point-like Yukawa particles, interacting via u(r) = Uoe~^ r / dm /(r /d m ), the 
state space is two-dimensional and spanned by Uq and A. In this space the melting line is a 
function Uq 1 (X) that has been determined by RKG. More precisely, they introduced an effective 
temperature T (kT in units of the Einstein phonon energy) which is related to Uq by T(A) = 
(2A 2 6'(A)/3C/ (A)/3) _1 where the function 6>(A) is given in Tab. I of 0. They then determined 
the melting line as f M (\) = 0.00246 + 0.000274A which is plotted as the solid line in Fig. (g].b). 
For colloidal systems, we recall that the Yukawa pair-potential of DLVO theory reads 



0u(r) 



ga 



L 1 + KeffCL 



2 



\b~ (2) 



with the effective colloidal charge Z e R and the effective screening parameter K e s . These effective 
(or renormalized) quantities, introduced to capture effects arising from the non-linearity of the 
PB equation jL8|, are both (unknown) functions of X, K c sa = fi(X) and Z e ff = f2(X). One 
way, among others, to approximately determine these functions is to use the PB cell model 



lq , |19| . Identifying the parameters in eq. (g) with Uq and A from the RKG calculation as 
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Figure 2: (a) Effective colloid-colloid pair-force as function of the distance between two colloids 
surrounded by 106 other colloids in a FCC configuration at rj = 0.03, Z = 1000, na = 0.2, 
As/a = 0.012. The dot-dashed line is the best fitting Yukawa interaction, (b) Z e R vs. Z, 
obtained from fitting Yukawa pair-forces to effective force curves in BCC (squares) and FCC 
(circles) configurations, for four salt concentrations as indicated and r\ = 0.03, Xb/cl = 0.012. 
The dot-dashed lines gives predictions of the PB cell model for comparison. 

0Uo = (Z c fie KcSCL / (1 + K c ga)) 2 XB/d m and A = K e gd m , one can then use the RKG phase-diagram 
in Fig. (|l|.b) to estimate whether a colloidal system, represented by a state point X, is solid or 
liquid. Note first that since the phase-behavior of a pure Yukawa system is determined just by 
Uq and A, the natural 3D state space of colloidal systems is reduced here to 2D. Secondly, it is 
clear that the data in Fig. (|].a) can be transfered to the phase-diagram of Fig. (|.b), only if the 
two functions fi(X) and f2(X) are provided. Then a data point at (Z,na) in Fig. (|l].a) can be 
transformed to (Z e s, K e ga) and subsequently plotted as (T,X) in Fig. (Oj.b). This requires two 
steps: (i) to calculate the effective pair- interaction (depending on the state X), which includes 
the original many-body interactions and, in general, is not a Yukawa potential and then try (ii) 
to fit these pair-interactions to Yukawa potentials. 

We proceed with the first task of finding an effective pair-interaction by performing step 
(i) to (iii) of the simulation algorithm described above for fixed colloidal configurations. The 
force Fab exerted by a particle A on another particle B is obtained as the difference between 
the total force acting on particle B with particle A present, F B , and the total force acting on 
particle B, F B , after removing particle A while leaving all other particles at their positions p0(| . 
Varying the position of particle A results in the effective force curve Fab{t) = F B {r) — F B . 
With this procedure the many-body interactions are folded into an effective pair interaction. If 
the true interactions in the system are pairwise additive, the resulting effective interaction is by 
construction identical to the true pairwise interaction potential, independent of the arrangement 
of the surrounding particles (all particles but A and B). This calculation has been carried out 
for FCC (N = 32 and 108) and BCC (N = 54) configurations at different salt concentrations 
between na = 0.1 and na = 2.0, bare colloidal charges between Z = 10 and Z = 4000, and 
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several volume fractions rj = 0.01, 0.02, 0.03. A typical effective force curve f(r) = a(3FAB{f) m 
a low-salt system (na = 0.2) is presented in Fig. (j2[a). Forces are multiplied by (r / a) 2 / (1 + K e gr) 
and plotted logarithmically so that the Yukawa interaction appears as a straight line with slope 
— K e ff. Obviously, for small distances the effective forces are Yukawa-like. However, around the 
mean colloid-colloid distance d m systematic deviations become visible which develop into a cut- 
off at r > dm. At the highest salt concentrations considered (na = 2.0), no such deviation occur 
in our calculation. A similar behavior has also been found in primitive model calculations ||. 

As a second step, the values of the effective parameters Z e Q and n e g are obtained by fitting 
the force curve to the derivative of eq. (||) at r < d m , that is, in the range where the force 
is Yukawa-like. Performing the fit for all calculated force curves yields the dependence of the 
effective parameters on the state of the system, i.e., the two functions fi(X) and f2(X). The 
dependence of Z c g on Z for different salt concentrations is shown in Fig. (§.b). The values 
obtained in FCC and BCC configuration are compared to the curves obtained from the PB cell 
model using the Alexander prescription |16L 19 1. The predictions of the cell model agree with our 
results in the high-salt case (na = 2.0), but not under low salt conditions {na = 0.2). Moreover, 
we find agreement between the FCC and BCC results at na = 2.0 but substantial differences 
at kcl = 0.2. As noted above, such a configuration dependence (ff^i-^) fi2 C {^)) °f the 
effective pair-interaction parameters is another fingerprint of the many-body interactions in the 
system next to the observed cut-off feature. A dependence of the pair-interaction on the colloidal 
configurations has also been found in [12|. In addition, attractive three-body forces have been 
explicitly shown to become effective in the low salt regime (na < 1) pi] . The range of the 
effective interaction in the high salt calculation (ko = 2.0), is about 1/K e ff ~ 0.5a, while for the 
low salt case {na = 0.2) it is roughly 2a. The mean distance d m in both cases is about 5a. The 
onset of the configuration dependence thus coincides with the size of the double layer becoming 
comparable to d m as pointed out in the introduction. 

Having determined the functions fi(X) and f2(X) the data for the melting line in Fig. (jl[a) 
can be transfered to Fig. (|].b), essentially by rescaling twice the x- and y-axis of Fig. (jl[a) 
((Z,Ka) — » (Z e ff,K e fja) — > (T, A)). Good agreement with RKG is obtained in our high-salt 
calculation na = 2.0, corresponding to large values of A, which is consistent with our finding 
that at high salt neither a cut-off behavior nor a configuration-dependence of the pair-potential 
could be observed. Obviously, in this salt regime the colloidal suspension can be represented by 
a Yukawa system quite well. However, reducing the amount of salt, i.e., decreasing A, we observe 
pronounced deviations from the RKG line occurring in the low salt regime. Again, this matches 
with the behavior of the calculated effective force curves at low salt, showing a configuration 
dependence but also the cut-off feature. Both observations suggest that the difference between 
the RKG melting line and ours is due to many-body effects. 

For practical applications it would of course be highly desirable to have a way of predicting 
the phase behavior of colloidal systems from a simple density-dependent pair-potential. The cut- 
off behavior observed in our effective force curves at low salt concentration and in the experiment 
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in H suggests to use a truncated Yukawa potential 

-Ar/d, 



u{r) = { ~" r ~ Tc (3) 

r > 

with a density -dependent cut-off r c oc d m to include the macro- ion shielding effect ||. With 
this model potential, we have carried out MD simulations and again determined the solid- 
liquid phase-boundary by the Lindemann criterion, computing the rms displacement for various 
combinations of Uq and A in FCC and BCC crystals. For r c = 3.07d m , the cut-off used for 



numerical reasons in [17], the RKG melting line is reproduced. Upon decreasing the cut-off, 
systematic deviations from the RKG melting line are observed occurring first at small values of 
A. The dotted line in Fig. (Jl|.b) is the melting line obtained choosing r c = 1.35<i m in a FCC 
configuration, while the dash-dotted line is the melting line in a BCC configuration choosing the 
cut-off radius r c = 1.50d m . These two cut-off radii correspond to the first neighbor shell in the 
FCC and to the second neighbor shell in the BCC configuration so that the number of neighbors 
that are included in the interaction is comparable in both cases. By comparing these results to 
our full PB-MD simulation it is seen that the density-dependent cut-off provides a considerable 
improvement over the plain Yukawa potential. Since the cut-off in this model pair-interaction is 
related in a simple way to the macro-ion shielding effect this further corroborates our conclusions 
about the importance of many-body interactions for the solid-liquid phase behavior of colloidal 
suspensions. 



Considering three-body forces for three colloids in an electrolyte, [21], one can show that 
the macroion shielding effect becomes less efficient with decreasing a. This implies that the 
observed difference between the RKG line and our line should become smaller with decreasing 
r\ = 47r(a/(i m ) 3 /3, which indeed we have observed in reference calculations at smaller volume 
fractions (r/ = 0.01 and r\ = 0.005). This is, first of all, a reminder of the fact, already discussed 
above, that the phase-diagram of colloidal suspension depends on three quantities, while the 
Yukawa phase-diagram is just a 2D representation of it. Secondly, this is the reason why the 



existing experimental data in [22, 23] can not be used to validate our results, as these experiments 
have been performed at even smaller volume fractions (10 -3 to 10 -4 ). For such low volume 
fractions we do not expect noticeable many-body effects. We would therefore suggest to perform 
experiments similar to the ones in [22, [23| at higher volume fractions comparable to r\ = 0.03 



considered in our simulation. 

A charge-stabilized colloidal suspension is often considered as a simple Yukawa liquid. This 
is not always correct. In conclusion, this paper shows that a Yukawa description of colloidal sus- 
pensions fails to predict the correct solid-liquid phase behavior under low-salt conditions where 
the range of the effective interaction is comparable to the mean distance and where many-body 
effects start playing a vital role. We have demonstrated this with the help of effective force cal- 
culations which have revealed a configuration dependence of the pair-interactions and deviations 
from Yukawa behavior, at r > d m and low salt. Model pair-potentials with density-dependent 
cut-offs have been shown to reproduce quite well the effects which many-body interactions have 
on the phase-behavior of colloidal suspension. We predict many body effects to increase with 
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increasing colloid volume fraction. 
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M. Brunner, Y. Chen, M. Brumen, D. Halozan and C Russ. 
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